#pragma OPENCL EXTENSION cl_amd_printf : enable
#define N_UNROLL  4  //unroll, future by mako


__kernel void step_v_self( 
    const int n_atoms,
    const float dt,
    __global float4 *f_e,      //force
    __global float4 *v_w       //velocity
    ) 
{
    int iglb = get_global_id(0);
    int sglb = get_global_size(0);
    float4 dt4 = (float4)(dt,dt,dt,0.f);
    int i = iglb;
    while(i<n_atoms){
      //for(int kk=0;kk<N_UNROLL;kk++){ //let compiler do unrolling 
        float4 vw = v_w[i];
        v_w[i] = vw + f_e[i]*dt4*vw.w; //w=1/mass
        i += sglb;
      //}end-unrolling
    }
}

__kernel void step_v( 
    const int n_atoms,
    const float dt,
    __global float4 *f_e,      //force
    __global float4 *v_w_in,   //velocity
    __global float4 *v_w_out  //velocity after dt
    ) 
{
    int iglb = get_global_id(0);
    int sglb = get_global_size(0);
    float4 dt4 = (float4)(dt,dt,dt,0.f);
    int i = iglb;
    while(i<n_atoms){
      //for(int kk=0;kk<N_UNROLL;kk++){ //let compiler do unrolling 
        float4 vw = v_w_in[i];
        v_w_out[i] = vw + f_e[i]*dt4*vw.w; //w=1/mass
        i += sglb;
      //}end-unrolling
    }
}


__kernel void step_r( 
    const int n_atoms,
    const float dt,
    __global float4 *v_w,      //veloctiy
    __global float4 *r_t_in,   //coordinates
    __global float4 *r_t_out   //coordinates after dt
    ) 
{
    int iglb = get_global_id(0);
    int sglb = get_global_size(0);
    float4 dt4 = (float4)(dt,dt,dt,0.f);
    int i = iglb;
    while(i<n_atoms){
      //for(int kk=0;kk<N_UNROLL;kk++){ //let compiler do unrolling 
        r_t_out[i] = r_t_in[i] + v_w[i]*dt4;
        i += sglb;
      //}end-unrolling
    }
}

__kernel void step_r_limited( 
    const int n_atoms,
    const float dt,
    __global float4 *v_w,      //veloctiy
    __global float4 *r_t_in,   //coordinates
    __global float4 *r_t_out,  //coordinates after dt
    const float4 max_dr       //max of allowed dr [mdx,mdy,mdz,0]
    ) 
{
    int iglb = get_global_id(0);
    int sglb = get_global_size(0);
    float4 dt4 = (float4)(dt,dt,dt,0.f);
    int i = iglb;
    while(i<n_atoms){
      //for(int kk=0;kk<N_UNROLL;kk++){ //let compiler do unrolling 
        float4 dr = v_w[i]*dt4;
        dr = fabs(dr)>max_dr ? max_dr : dr;
        r_t_out[i] = r_t_in[i] + dr;
        i += sglb;
      //}end-unrolling
    }
}

